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ABSTRACT 

We present numerical simulations of penetrative convection and gravity wave excita- 
tion in the Sun. Gravity waves are self-consistently generated by a convective zone 
overlying a radiative interior. We produce power spectra for gravity waves in the ra- 
diative region as well as estimates for the energy flux of gravity waves below the con- 
vection zone. We calculate a peak energy flux in waves below the convection zone to be 
three orders of magnitude smaller than previous estimates for m=l. The simulations 
show that the linear dispersion relation is a good approximation only deep below the 
convective-radiative boundary. Both low frequency propagating gravity waves as well 
as higher frequency standing modes are generated; although we find that convection 
does not continually drive the standing g-mode frequencies. 

Key words: Sun:mixing, internal gravity waves, convection 



1 INTRODUCTION 

Convection bounded by a stable region is important in all 
aspects of stellar evolution. The depth to which convec- 
tive motions can overshoot into an adjacent stable layer 
has important consequences for many stages of stellar evo- 
lution, from dredge up in lat e stages of stellar evolution 
iKippenhan a nd Wcigert"l994') to nucleosynthesis in C las- 
sical Novae ijWooslev 1986 : Kercok c t al 199& : .Kercek et alJ 
119991: lAlexakis "et "al.1i2004i) . In the Sun, the overshoot region 
at least partially overlaps the tachocline (the shear layer at 
the base of the convection zone where the rotation changes 
from being differential in the convection zone to solid body 
in the radiation zone and the likely site of toroidal field gen- 
eration in the dynamo process). Therefore, understanding 
the extent of the overshoot is necessary for understanding 
the nature of the tachocline; is it a violent region, constantly 
bombarded by overshooting plumes, or is it a laminar region 
just below the overshoot? The answer to this question can 
constrain models for the solid body rotation of the solar 
interior. 

In addition to the uncertainty of penetrative convec- 
tion, which has plagued stellar evolution theory, there are 
other noteworthy aspects of convective-radiative boundaries. 
Convection adjacent to a radiative region generates inter- 
nal gravity waves which can have dynamical affects. In the 
Earth's atmosphere, the nonlinear interaction between con- 
vection, shear and gravity waves is the cause o f the quasi- 
biennial oscillation (QBO) iBaldwin et aljEooj) : this shear- 
gravity wave interaction has been demonstrated in a remark- 



E-mail:trogers@pmc. ucsc.edu 



able experiment bv lPlumb and McEwanI il978tl . The same 
process responsible for th e QBO has been invoked to ex plain 
the solar magnetic cycle iKumar. Talon fc Zaiinlligggll . 

Internal gravity waves can share angular momentum 
with the mean fiow when they are attenuated. The trans- 
port and deposition of angular momentum by gravity waves 
is one explanation for the solid body rotation observed in 
the radiative interior of the Sun (iKumar fc Quataerdll99'^ : 
IZahn. Talon, fc Matialligg?!) . However, a vast literature in 
atmospheric science tells us that gravity waves have an 
anti-diffusive nature; they enhance local shear rather than 
smooth it. In addition, internal gravity waves have also been 
implicated in the mixing of species (|Presjll98lD. and in 
order to e xplain the lithium depletion dSchatzmanl Il993l : 
ICarcia-Lopc z and Spruit 1991j). 

Most studies of convection bounded by a stable 
layer have been focused on the penetration depth. The 
depth that convection overshoots into an adjacent ra- 
diative region has been studied extensively by numerical 



simulations in two I Massaeuer et al.||l984^ [Hurl burt et alJ 
19861: iHurlburt et allKbg^ iRoeers fc GlatzmaieJl2005^ and 



three dimensionl llBnnnmeU^^l .l l2002l:ISaikia et alJl2000H . 

In general, these simulations find that the penetration 
depth is a sensitive function of the stiffness of the sub- 
adiabatic layer and find that it is also sens itive to 
the Peclet number (^) jBrummell et alJ l2002h . Early 
2D studi es with "s oft" interfaces (mild subadiabaticity) 
iHur lburt et al.llT98&: Hurlburt ctHligg^ found that over- 
shooting convection could ext end the adiabatic re gion. 
Three-dimensional simulations ijBrummell et alJl2002^ and 
simulati ons with very stiff (large su badiabaticity) radiative 
regions llRoeers fc GlatzmaieJl2005l) find no extended adia- 
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batic region, in agreement with helioseismic observations. 
However, in all of the previous studies a constant (with 
height) subadiabaticity was used and none were as stiff as 
the Sun. Analytic models generally predict extended adia- 
batic regions, but rely on highly uncertain parametrizations. 
The penetration depth in the Sun and its dependencies will 
be discussed in greater detail in a forthcoming paper. 

While internal gravity waves have b een studied ana- 
lytically for some time, e.g. IPressI lll98 J) , they have not 
received much attention in numerical simulations of stel- 
lar interiors. In the few instances where they have been 
simulated they have not been analyzed in detail and have 
been unlike those in the Sun because the subadiabaticity 
was much less than solar. However, iKiraga et alJ (l2003f) 
did extensive numer ical studies in 2D, similar to those by 
iHurlburt et al l)l994) and analyzed the wave spectrum pro- 
duced and compare d it to previous analytic predict ions for 
g-mode spectra bv iGarcia-Lopez and Spruid il99ll) . They 
found that the spectrum of waves produced in numerical 
simulations was more broadband in frequency and that nu- 
merical simulations always predicted larger fluxes in waves 
than analytic studies. The overestimate of energy in waves 
in that study was probably due to the mild subadiabatic- 
ity used in those calculations. The subadiabaticity of the 
stable re gion is extremely important in these types of sim- 
ulations iRogers fc Glatzmaier! |2'005^ : the mild values used 
in all previous simulations allow plumes to penetrate much 
more deeply, causing larger nonlinearities and velocities than 
would otherwise occur. 

However, since internal gravity waves may be main- 
taining the solid body rot ation of the solar interior 
l)Zahn. Talon, fc Matiaslll997^j_ may be driving a QBO l ike 

oscillation (coined SLO bjHSIs^^^^ii^SSsil-llQS^LJ. 

at th e base of the convection zone " |KumM^^kmfc^aSi] 
[l999|) and may someday be observed with helioseismology 
iA^pourchaux et alj|200ol) . we have begun a detailed study 
via numerical simulations of their excitation and evolution. 
The way in which gravity waves redistribute angular mo- 
mentum and the likelihood of their observation depends on 
the spectrum of gravity waves produced by the overlying 
convection, their nonlinear interaction with convection and 
rotation, and on their radiative damping. Here we discuss 
the spectrum of gravity waves produced by overlying con- 
vection and penetration, radiative damping and the energy 
flux of waves beneath the convection zone. The impact grav- 
ity waves have on the angular velocity in the radiative in- 
terior will be discussed in a forthcoming paper. We discuss 
the numerical method in section 2, convection in section 3 
and gravity waves in section 4. 



2 NUMERICAL MODEL 
2.1 Equations 

We solve the Navier Stok es equations w ith rotation in the 
anelastic approximation l|Goughl Il969l) for an ideal gas. 
The anelastic approximation is appropriate when the fluid 
flow is sufficiently subsonic and the perturbations relative 
to the mean (reference) thermodynamic state are small. 
The radially (r) dependent reference state is taken from a 
ID solar model (Christensen-Dalsgaard private communica- 



tion). Unlike previous formulations of the anelastic equa- 
tions jGlat zmaie il ll985l: iRogers et alJliool) . here our refer- 
ence state temperature gradient is not adiabatic, and we 
choose to use the temperature perturbation as our work- 
ing thermodynamic variable, instead of the specific entropy. 
This formulation has advantages an d disadvantages com- 
pared to the standard formulation iRogers fc GlatzmaieJ 
l2005f) . The equations are solved in 2D cylindrical coordi- 
nates, with the radiation zone extending from O.OOIRq to 
O.7I8R0 and the convection zone occupying the region from 
O.7I8R0 to 0.9Rq. This geometry represents an equatorial 
slice of the Sun that includes most of the radiative and con- 
vective regions. For numerical reasons we put another stable 
region from 0.9Rq to O.93R0. The exclusion of ^% of Rq 
around the centre has little effect since gravity waves inter- 
nally reflect where their frequency equals the Brunt- Vaisala 
frequency, which vanishes at the centre. 

We solve the following anelastic equations for perturba- 
tions to the reference state: 



V • = 0. 



dt 



+ {v ■ V)v = -VP - Cgf + 2{v X !^)- 



(1) 



(2) 



dT dT — 

^ + (^T. V)T = -vr{^ - (7 - l)Th,) + 

(7 - l)ThpVr + 7K[V^r + [hp + '^«)^] + 

UT Cy 



(3) 



Equation (1) represents the continuity equation in the 
anelastic approximation, where p is the reference state den- 
sity and V is the fluid velocity. Equation (2) is the momen- 
tum equation, assuming a constant dynamic viscosity, pv. 
The gravitational acceleration, -"gr is taken from the ID so- 
lar model. The rotation rate, O is set to the mean rotation 
rate of the Sun 2.6 x 10~^ ra d/s. The reduced pressure, P 
iSraginskv and Robert ^119951) is defined as 



+ U 



(4) 



where p is the pressure perturbation and U is the gravita- 
tional potential perturbation. The co-density perturbation 
is defined as: 

C = -i(T+4^p). (5) 

Equation (3) is derived from 

dT_ _ dT dp dT _ 

dt dp'" dt'^^ds'" dt ' 

where here the thermodynamic variables are the full vari- 
ables (sum of the reference state, denoted by an overbar, 
and the perturbation) and an ideal gas is assumed. The 
thermal expansion coefflcient a = =, the adiabatic expo- 
nent 7 = ^ and K, the thermal diffusivity are all functions 
of radius. The inverse density scale height and thermal dif- 
fusivity scale height are defined as: 



V-v-i V-(cp/9KVT')(6) 

a pCv 
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h -^^ 

p ar 



and 



, 1 dK 

hK = --r 

K ar 



(7) 



The first term on the right hand side of Equation (3) 
can also be written as 



^.dlnT dlnT 



(8) 



which illustrates how this term depends on the local diflFer- 
ence between the nonadiabatic reference state and an adia- 
batic state. The last two terms in Equation (3), which are 
only a function of the reference state, are assumed zero at 
the beginning of our simulations. We neglect viscous heat- 
ing for these low viscosity simulations. The top and bottom 
boundary conditions are isothermal, stress-free and imper- 
meable. 



2.2 Numerical Technique 

The technique use d her e is similar to that of 
iRogers fc Glatzmaieil (|200^, but for a rotating disk 
instead of a box. For numerical simplicity we take the curl 
of the momentum equation, leaving us with a vorticity 
equation: 

duj V, Q dT 1 dT dp _„2 /„x 

-+(«.V). = (2n-f.)^«.-^---=--^ W.(9) 

where the vorticity, iu — V x v is in the z direction. 
We have neglected viscous terms involving ^ . We calculate 
the pressure term using the longitudinal component of the 
momentum equation (2): 

1 dp dve - „ 2^. hpdvr 

where ve and Vr represent the longitudinal and radial ve- 
locity components, respectively. The time derivative of ve is 
from the previous timestep and we neglect the perturbation 
to the gravitational potential in (10). 

These equations are solved using a Fourier spectral 
transform in the longitudinal (9) direction and a finite dif- 
ference scheme on a non-uniform grid in the radial (r) di- 
rection. Time advancing is done using the explicit Adams- 
Bashforth method for the nonlinear terms and an implicit 
Crank-Nicolson scheme for the linear terms. 

We do a 21st order polynomial fit to the reference state 
density (p), temperature (T), mass (which gives us gravity, 
g), opacity (k) and 7 given in the ID solar model. Figure 1 
shows density and temperature as a function of radius in our 
model (solid line) compared to the standard model (dotted 
line). Analytic derivatives of these (fitted) values give us 
inverse scale heights {hp,hT and /i^) and the temperature 
gradient. 

The subadiabaticity in the radiative region is then de- 
fined as: 



dT 
dr 



{1-1)KT 



(11) 



where the second term in equation (11) is the adiabatic tem- 
perature gradient. Figure 2 shows the effective Brunt- Vaisala 
frequency (subadiabaticity), defined in (16) as a function of 
radius. The superadiabaticity in the convection zone is spec- 
ified to be a constant, 10~^. The thermal diffusivity is that 



given by the solar model, multiplied by a constant factor for 
numerical stability: 

-3 



K = kapmult ■■ 



16ar 



(12) 



where a is the Stefan-Boltzman constant and the multiply- 
ing factor, kapmult, is 10^ (and therefore, the convective 
flux is 10^ larger than Solar). This gives us the proper ra- 
dial proflle of the radiative diffusivity, albeit increased by 
a large factor for numerical stability; this is a "turbulent" 
diffusivity. The viscous diffusivity is specified to be: 



const 



(13) 



keeping the dynamic viscosity constant. Since both V and 
K vary with height, the Pr (Prandtl number = =) varies 
with height from 10^^^ (at base of the convection zone) 
to 0.7tcz (at top of convection zone) and is 10~^ near the 
centre. Therefore, the Ra (Rayleigh number = ^^^=-2-, 
where D is the depth of the convection zone and AVT is 
the superadiabatic temperature gradient) also varies from 
~ 10* at the base of the convection zone to ~ 10^ at the 
top of the convection zone. The Ekman number (Ek= 2(10'^ ) 
varies from W^q^ to 10^^^. The resolution used in these 
models is 2048 longitudinal zones x 1500 radial zones, with 
500 radial zones dedicated to the radiative region. 



3 CONVECTION BOUNDED BY A STABLE 
LAYER 

Here we discuss the convection in this model. In order to 
illustrate the differences between convection bounded by a 
radiative region and convection bounded by a hard bound- 
ary we also present a model that has an impenetrable bot- 
tom boundary. The model with an impenetrable boundary 
has the same reference state as the model with a stable re- 
gion and therefore, has the same Ra, Pr and Ek. We note 
here that the simulations with impenetrable boundaries usu- 
ally require greater spatial resolution near the boundaries 
and smaller timesteps because the hard boundaries deflect 
fast descending (ascending) plumes into shallow, high-speed 
horizontal flows. This problem is more severe at the top 
of the computational domain because the density is lower 
there, resulting in larger velocities. This is one reason we 
have added a very shallow stably stratified region at the top 
boundary. Another is that it is more realistic to have rising 
plumes stopped by a stable layer rather than an artificial 
hard boundary. 

Snapshots of the vorticity in the two models are shown 
in Figure 3. There are many differences immediately obvious 
in this figure. Most noticeable is the nature of the convec- 
tive cells. In the model with an impenetrable (hard) bottom 
boundary, there are fewer convective cells and those cells are 
more orderly and appear more laminar. Because descending 
and ascending plumes are diverted horizontally by the hard, 
stress-free boundary nearly all of their vertical kinetic en- 
ergy is converted into horizontal kinetic energy. These (rel- 
atively) large horizontal velocities allow the flow to travel 
further before cooling (heating) enough to fall (rise). In the 
model bounded below by a stable layer, a small fraction of 
the vertical energy in descending plumes is transferred to 
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the stable region generating gravity waves. Because of this 
transfer of energy and less rigid horizontal trajectories, hor- 
izontally diverted fluid does not travel as far before cooling 
(heating) and sinking (rising). The time dependent convec- 
tive penetration in the model bounded by a stable region 
also results in more chaotic interaction between convective 
cells, albeit with lower overall kinetic energy. 

Figure 4 shows the kinetic energy density (wg ■^- Vr), 
horizontal kinetic energy flux vg * (i;| + v^) and vertical ki- 
netic energy flux Vr * (wg -\- v'^), all scaled to the maximum 
values in the impenetrable case. It is obvious in this flgure 
that in the convective region there is significantly more en- 
ergy (by 2-3 orders of magnitude) in the model with hard 
boundaries. This is due almost entirely to a difference in 
magnitude of horizontal velocity (4b, note magnitudes), as 
explained above. 

In models bounded by a stable region a fraction of the 
kinetic energy is shared with the radiative region. Descend- 
ing plumes overshoot the convective-radiative boundary and 
energy is transferred to gravity waves. The simplest estimate 
of the penetration depth is given where the kinetic energy 
density drops to 1% of its maximum value in the convection 
zone. In this model this results in a penetration depth of 
approximately 0.1 pressure scale heights, but no extended 
adiabatic region. Several things can affect this penetration 
depth. First, the "stiffness" of the subadiabatic layer has 
been shown to be of utmost importance in several previ- 
ous numerical studies |ll urlburt et alJlQSStlBrummell et alJ 
l2002l: iRoeers fc GlatzmaieJ l2005l) . As stated above, this 
model has the (realistic) subadiabatic structure of the Sun 
and hence, a very high "stiffness" , which makes the penetra- 
tion depth small. The dimensionality of the model also has 
an effect; 3 D models have smaller penetration depths than 
2D models llBrummell et al.ll2002h . Therefore our estimate 
may be an upper limit. In addition, rotation has been shown 
to hi nder penetration iBrummell et alJ l2002l : lJulien et alJ 
Il996t ), while the effect of the Ra yleigh number is still un- 
clear jRogers fc GlatzmaieJf2005ll and may depend on the 
dimension. 

Kinetic energy spectra in the convection zones of the 
two models also differ. Figure 5 shows the kinetic energy as 
a function of horizontal wavenumber at a radius half way 
through the convection zone in both models. As illustrated 
in Figure 4, and now again in 5, the overall kinetic energy in 
the convection zone is larger for the model with an impene- 
trable lower boundary (shown as the dotted line in Figure 5), 
especially at low wavemode (large scales). Both models are 
fit rather well with a power law. While fc~"^ spectrum is 
expected in 2D turbulence, the geometry of the fiow has not 
been previously considered, nor the effect of adjacent stable 
regions (both models have stable regions on top) . Fully con- 
vective simulations of convective motion in a 2D disk give 
a variety of scaling laws (between k^^ and fc~^) depend- 
ing on the Ekman number, the density stratification and 
the depth at which the spectrum is calculated. The model 
with hard boundaries shows less of a dissipation range and 
some build-up at small scales. The model bounded by a sta- 
ble region shows a dissipation range and no such build-up 
at large wavemodes. This difference is probably due to the 
issues mentioned above: vorticity generation and large hor- 
izontal velocities along the (impenetrable) boundary. 

In summary, both the energetics and the character of 



convection are changed when a bounding stable layer is used 
instead of the traditional (artificial) impenetrable boundary. 



4 GRAVITY WAVES 

Gravity waves can be generated by Reynolds stresses in a 
convective region overlying a radiative region as well as by 
descending plumes that overshoot the convective-radiative 
boundary. The gravity waves generated by such processes 
have been studied previously in analytic studies. In partic- 
ular, iGoldreich. Murrav fc KumaJ (Il994l) studied the spec- 
trum of gravity waves generated b y Reynolds stresses in the 
overly ing convective region and iGarcia-Lopez and Spruiti 
il99ir) studied the spectrum of gravity waves generated by 
overshooting plumes. Both models are analytic and assume 
some combination of mixing length theory (MLT) and a Kol- 
mogorov spectrum for turbulence. 

Gravity waves have the ability to mix species and trans- 
port angular momentum and thus, may have profound ef- 
fects on the dynamics of the radiative zone in the Sun. In the 
following, we discuss results from both linear and nonlinear 
numerical simulations. First, however we briefiy review basic 
properties of gravity waves in the asymptotic r egime with- 
out ro t at ion (for a more thorou g h revi e w see, e. g.lUnno et all 
^1989^ : lchristensen-Dalsgaardl ll2002ll : iGouehl lll99ll) ). The 
asymptotic, linear dispersion relation for pure gravity waves 



Nkh 
io = — - — 
k 

where kh represents the horizontal wavenumber and 

k = {kl + kl)^ 



(14) 



(15) 



where kr is the radial wavenumber. The Brunt- Vaisala fre- 
quency, N, represents the upper limit on the frequency of 
gravity waves. It depends on the thermal stratification of 
the fluid and is defined as: 



T dr 



{l~l)hpT) 



(16) 



From the basic linear dispersion relation (14) it is seen that 
the frequencies of gravity waves depend only on the gas 
stratification and the direction of phase propagation. From 
the above dispersion relation it can also be shown that the 
group velocity and phase velocity are perpendicular, a pe- 
culiar trait of internal gravity waves. 

In addition to the standing g-modes (in the frequency 
range of 100/.iHz to 500/iHz), propagating waves from the 
overlying convection and penetration are produced. These 
propagating gravity waves, generally with much lower fre- 
quencies, have been suggested as the source of the oscil- 
lating shear layer at the base of the convection zone, the 
net angular momentum transport in the radiative interior 
CZahn. Talon, fc Matias"l997l: iKumar. Talon fc Zahnlll99gl : 
Talon. Kuinar,_fc Zahn 20 02i) a nd the mixing of species 
IIGarcia-Lopez and Spruiti Il99ll) . The spectrum of these 
waves has been predicted in those analytic studies. However, 
here we present the first self-consistent nonlinear simulation 
and analysis of these waves in the Sun's radiative interior. 
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4.1 Linear Gravity Waves 

While our main simulation here is nonlinear in both the 
convective and radiative regions of the computational do- 
main, we have computed linear cases for comparison with 
previous asymptotic results and to illustrate the differences 
between linear and nonlinear solutions. We first discuss the 
linear calculations. Figure 6 shows a snapshot of the vor- 
ticity in the radiation zone produced with a linear version 
of our model, excluding the convection zone and initially 
perturbed at several horizontal wavenumbers. While there 
are several different horizontal wavenumbers depicted in this 
image, their modes do not interact. Thi s is similar to that 
predicted from ray-tracing llGoughlll99ll) . 

In our partially linear simulation, the full set of non- 
linear equations is solved in the convection zone, but the 
nonlinear terms in equations (3), (9) and (10) are decreased 
linearly between O.718i?0 (the radius of the transition be- 
tween convective and radiative layers) and O.69-R0. There- 
fore, below O.69-R0 only the linear set of equations is solved. 
A snapshot of vorticity looks similar to that in 8a (the fully 
nonlinear case). However, there is a remarkable difference 
between this gravity wave pattern (seen in 8b) and that 
in figure 6. When constantly being driven by penetrative 
convection at many different frequencies and wavenumbers 
g-modes form a much simpler spiral pattern in the radiative 
interior. 

A simple reason for this behaviour can be explained as 
followed. If convection is considered as nearly equally spaced 
convective cells, descending plumes drive similar wavenum- 
bers, albeit slightly out of phase. The nonlinear wave-wave 
interactions generate new waves with wavenumbers equal 
to both the sum and difference of the original waves. The 
high wavenumber waves are damped in a very short dis- 
tance; wherease the small wavenumber (long wavelength) 
waves survive through most of the radiative interior. This 
wave- wave interaction only occurs in nonlinear simulations. 
Figure 6 has no nonlinear interactions and thus all initially 
excited wavenumbers remain. 

Frequencies of gravity waves in the solar interior in the 
linear, non-dissipa tive, non-rotating case have been calcu- 
lated previously ( llUnno et allll98'9riChristensen-Dalsgaardl 
|20o3) and generally range from lOO/xHz at low horizontal 
mode (large wavelength) to 500/iHz for higher horizontal 
mode (small wavelength) depending on the radial order. Fig- 
ure 7 shows the power spectra of gravity waves in the radia- 
tive interior in our partially linear simulation shortly after 
the simulation was started. Significant power occurs at the 
base of the convection zone, up to the highest frequencies for 
the low modes and up to about 200pHz for the high modes 
(figure 7a). At this radius the fluid motion is fully nonlin- 
ear and the expected broadband nature of the convection is 
seen. Note that at the convective-radiative boundary there 
is power in frequencies much higher than the typical con- 
vective turnover frequency of 2-lO^Hz because descending 
plumes pummel the stable region at many locations slightly 
out of phase. Power is not concentrated at the convective 
turnover frequency as is often assumed. 

Moving into the stable layer, figure 7b shows the dis- 
persion relation for gravity waves at 0.575 Rq. The peak 
power in figures 7b and c is three orders of magnitude 
smaller than the peak power shown in figure 7a. There is 



(significant) power at these smaller radii only out to hor- 
izontal mode 25, indicating the rapid decrease in energy 
due to radiative damping. In figure 7b there is power in 
both low frequency propagating g-modes generated by the 
overlying convection and penetration (below lOO/xHz) and 
high frequency (above lOOpHz) standing waves. Ridges sim- 
ilar to those calculated by asymptotic analysis (Christensen- 
Dalsgaard 1986) for high radial order are apparent. At any 
given height there are several radial orders contributing to 
the power at that horizontal mode and frequency, giving rise 
to the width of the ridges seen in this dispersion relation. At 
radius 0.360 Rq (figure 7c), it is seen that power at higher 
modes is diminished relative to figure 7b and that the lowest 
frequency /high wavenumbers have been completely damped 
(see section 4.2.2). Again, there is significant energy in the 
propagating waves at frequencies much higher than the con- 
vective turnover time. 

4.2 Nonlinear Gravity Waves 

Figure 8a shows a snapshot of the vorticity produced with 
our fully nonlinear model for the full computational domain 
(as in figure 3b). Figure 8b shows the temperature pertur- 
bation just in the radiation zone to illustrate the salient fea- 
tures of the gravity waves produced; waves spiral toward the 
centre due to rotation and geometry. While both small and 
large horizontal modes are produced at the interface only 
small mode numbers survive at the centre with substantial 
amplitude. 

Figure 9 shows the power spectra for our fully nonlinear 
case after the convection has reached a statistically steady 
state. The peak amplitude in figure 9 is scaled the same 
as in figure 7. There is little energy in the high frequency 
standing waves. That is, there is little power at frequencies 
above lOO/iHz for modes greater than 10, in figures 9b and c. 
The lack of standing modes after the convection has reached 
steady state is seen in both linear and nonlinear calculations, 
indicating that the energy in these frequencies is radiatively 
damped and is not continually generated by the overlying 
convection. 

The appearance of the high frequency standing waves 
is a numerical artifact and we show them here purely as a 
proof that the model can reproduce the asymptotic limit 
(given a large enough perturbation). The simulations are 
started with an initial perturbation to the convection zone. 
This initial perturbation has much higher amplitude than 
steady state convection. The large amplitude perturbation 
hits the stable region hard enough to excite all resonant 
frequencies. Once the convection has reached a statistically 
"steady" state, high frequencies still exist in the convection, 
but at a lower amplitude. These lower amplitudes are unable 
to drive the high frequency waves. 

The (dis)appearance of the high frequency standing g- 
modes was investigated in several different models. We ran 
a series of models including: partially linear, fully nonlin- 
ear, increased rotation rate, decreased rotation rate, larger 
and smaller diffusion coefficients. In general any "large" per- 
turbation to the momentum equation, like increasing or de- 
creasing the rotation rate by a large value, or turning off 
nonlinear terms gives rise to the high frequency standing 
modes; however, these modes disappear on a timescale of 
^2-5xl0^s. It appears that statistically "steady" convec- 
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tion does not provide a large enough perturbation to con- 
tinually generate these modes. 

In the Sun the thermal diffusivity is 10^ times smaller 
than in our simulation and therefore, one would expect these 
high frequency waves to damp on the order of 10^^ s rather 
than the 10^ s seen in our simulation. However, this is still 
a relatively short time; and therefore these high frequency 
standing waves likely do not exist in the Sun. In addition, 
it has been found in previous studies that two dimensional 
simulations produce larger perturbations to the underlying 
stable region than 3D simulations. In 3D these modes would 
be excited with lower amplitude and hence, would be even 
less likely to exist. In order for these modes to exist, they 
will have to be continually generated by something other 
than steady convection. 



10 and 20/iHz,m=l waves dropping by only around 3 or- 
ders of magnitude. Also as expected, higher values of hor- 
izontal mode numbers are damped more than m=l and 
their energy in the lower frequencies is damped over a 
shorter distance. For comparison, the dot-dashed line shows 
the analytic prediction for the damping of the m=l wave 
jKumar. Talon fc Zahn| [l9991. Using amplitudes for the ki- 
netic energy at the convective-radiative interface from these 
simulations, w e damp them with depth usi ng the analytic 
prescription in lKumar. Talon fc Zahia |^9^; there is a sub- 
stantial difference, particularly for the lowest frequencies. In 
figure 11 it is seen that (1) only the very low m waves make 
it to the centre with any appreciable energy, (2) all three 
frequencies are excited with remarkably similar amplitudes 
and (3) m=l looks like a standing wave. 



4-2.1 Radiative Damping 

As waves propagate in a diffusive medium they are ra- 
diatively damped. This damping leads to the transfer of 
their angular momentum, positive or negative, to the dif- 
ferential rotation. The radiative damping of waves depends 
on the frequency and wavenumber of the waves. Low fre- 
quency and high wavenumber (small wavelength) waves are 
damped more quickly (over a shorter distance) than high 
frequency/low wavenumber waves. Because thermal diffu- 
sion has more time and larger temperature gradients for the 
former. 

In this geometry group velocity is directed inward and 
phase velocity is outward. The inward propagation of wave 
packets, indicative of group velocity, can be seen in figure 10. 
This figure shows the temperature perturbation as a func- 
tion of radius for a given horizontal mode (in this case m=20 
and 100, where the horizontal mode number m is defined as 
fc^r) for two sequential times, with the dotted line represent- 
ing the earlier time and the solid line representing a later 
time (lO'^s later). It can be seen that wave packets move 
inward in time. As wave packets move inward they are ra- 
diatively damped and therefore, their amplitudes decrease. 
In figure 10, not only does the lower m case have higher 
amplitude overall, it is damped less in time relative to the 
higher m case. In animations of these data a group velocity 
can be (crudely) measured and we find that at these radii 
(far away from the transition) this group velocity agrees well 
with that predicted by the linear dispersion relation. While 
we are able to pick out wave packets and a group velocity for 
higher values of the horizontal mode m, we see only standing 
waves in lower values of m. It is likely that there is a prop- 
agating wave component, however, it has significantly less 
amplitude than the standing waves and is therefore, hard 
to discern. This has implications for energy in propagating 
waves (section 4.2.2). 

The radiative damping of waves can also be seen in 
figure 11. In that figure we show the kinetic energy den- 
sity in waves as a function of radius for three different fre- 
quencies and horizontal mode numbers. The energy at the 
convective-radiative interface varies by no more than an or- 
der of magnitude across the three different frequencies and 
modes shown. However, moving radially inward it is seen, 
as expected, that the lowest frequency waves are damped 
more rapidly with depth, with the 2/iHz, m=l wave am- 
plitude falling by nearly 6 orders of magnitude and the 



4-2.2 Energy in g-modes 

Gravity wave energy and fiux is of great interest for mixing 
and angular momentum transport in the solar radiative re- 
gion. The time averaged kinetic energy density as a function 
of radius in the stable region is shown in figure 12. The large 
kinetic energy just below the convection zone is due to con- 
vective overshoot. The local peak in energy density at the 
centre may be due to the focusing affect discussed by Press 
(1981). We can also calculate the wave kinetic energy as a 
function of horizontal wavenumber (fc^) and frequency (/) 
at a particular radius: 

E{iU, kh) = ^p(ve(Lo, khf + Vr{uJ, khf). (17) 

This can be compared with the energy calculated using (14) 
and (15) : 

E{lJ, kh) = ^-p{-fvr[uj, kuY (18) 

where kh ~ —. The (dis)agreement of these two equations 
gives us an estimate of the applicability of the linear disper- 
sion relation at different radii, wavenumbers and frequen- 
cies. Figure 13 shows these two relations at two radii and 
two horizontal modes as a function of frequency. It is seen 
in that figure that the two formulations do not agree near 
the convective-radiative interface, with the linear disper- 
sion relation underestimating the energy by 2-1- orders of 
magnitude for both m=l and m=10. Well away from the 
convective-radiative boundary, the linear dispersion relation 
is a good approximation, particularly at m=10 where the 
two formulations match exactly (the dotted line is under the 
solid line). It is expected that the linear dispersion relation 
would break down near the convective-radiative interface, 
given that overshooting plumes cause significant nonlinear- 
ity. 

Since the linear dispersion relation does not hold near 
the convective-radiative boundary it is difficult to justify 
the use of the group velocity calculated from that disper- 
sion relation at those radii. Therefore, in order to estimate 
the energy ffux in waves at the base of the convection zone 
(defined as energy density in waves times the group veloc- 
ity times the surface area), we use the group velocity at a 
deeper radius, where the linear dispersion relation is a good 
approximation. 
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The wave flux calculated this way is shown in figure 
14. The peak flux here (m=l,low frequency) is three or- 
ders of magnitude smaller than the peak wave flux pre- 
dicted in lKumar. Talon &i Zahn (1999) for 1=1. However, in 
our simulations this flux is distributed evenly in frequency 
for ah (relatively low) wavemodes unlike the calculations 
by iKumar. Talon fc Za hn ( 1999) where the flux drops ofl' 
rapidly with frequency. This broadband distribution of en- 
ergy produced is seen in figure 1 1 and was found in the 
simulations bv lKiraea et al.l ll2003h . 

These estimates should be considered rough at best, as 
the group velocity at radii close to the convective-radiative 
boundary is very uncertain. The kinetic energy of the fluid, 
calculated using (17) is likely not just wave energy. How- 
ever, this energy, while not completely wave energy, still 
contributes to the angular momentum transport and mix- 
ing below the convection zone. In addition, these estimates 
should be considered as upper limits on the wave flux for 
a couple of other reasons; (1) in two dimensions energy is 
transferred to the stable region more readily than in 3D, be- 
cause of the alignment of convective cells and the "flywheel" 
motion that tends to occur in 2D convection because of the 
inverse cascade, (2) our thermal diffusivity is signficantly 
larger than the solar value (by a factor of 10^), causing the 
convective velocity to be larger than those expected in the 
Sun and therefore, the wave flux could also be larger and (3) 
most of the energy in low frequency/ low mode waves are in 
the standing wave component rather than the propagating 
wave (figure 11). 



5 CONCLUSIONS 

We have presented self consistent calculations of gravity 
wave excitation by overlying convection and penetration. 
These simulations produce power spectra of the gravity 
waves as a function of radius. We find that both high fre- 
quency standing g-modes and low frequency propagating 
modes are excited depending on the amplitude of the per- 
turbation. In general the standing waves are not continually 
driven by the convection and only larger perturbations to 
the momentum equation produce these modes. 

Nonlinear effects need to be considered in the radia- 
tion zone for several reasons. These effects broaden the fre- 
quency ridges in the dispersion relation, they allow for trans- 
fer of energy from high frequency to low frequency modes. 
In addition, just beneath the convection zone, where plumes 
overshoot, nonlinear interactions increase the energy there 
by more than two orders of magnitude over what the lin- 
ear dispersion relation would predict for energy in waves. 
While this energy may not be completely wave energy it 
still contributes to the dynamics of the region just below 
the convection zone. 

The low frequency modes are driven by the overlying 
convection and penetration. Typical frequencies of these 
waves can be ten times larger than the convective turnover 
frequency, which varies from 2-lOpiHz. We also find that the 
energy in these waves is distributed rather evenly in fre- 
quency. The steep drop off in energy at larger frequencies is 
not observed. 

While these are self-consistent calculations of gravity 
waves in a realistic solar simulation, there are many short- 



comings that still need to be addressed. First is the three 
dimensional effect on the waves and probably more impor- 
tantly on the nature of the convection which drives the 
waves. In addition, any numerical simulation will be more 
laminar than the Sun. Presumably, a more turbulent convec- 
tion zone should provide a broader spectrum (in frequency 
and wavenumber) of waves. What effect this will have on 
the spectrum and dissipation of the waves is unclear. 
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Figure 1. Reference state density and temperature, taken from 
the ID solar model. Both fitted values (solid line) and a<;tual 
values (dotted lines) are shown (dotted line is directly under the 
solid line). 
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Figure 2. Brunt- Vaisala frequency squared as a function of radius 
in the stable region. Calculated using equation 16. 
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Figure 3. Vorticity (in ^-^) for a convection model with im- 
permeable boundaries (a) and for convection model with a sta- 
ble layer below (b), blue represents negative vorticity while 
white/yellow represent positive vorticity. 
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Figure 4. Kinetic Energy Density (2^) (a) & (d), Horizontal 
Kinetic Energy Flux (b) & (erand Vertical Kinetic En- 

ergy Flux ( ) (c) & (f)- AH plots are scaled to the pealj value 
in the impermeable case. The qualitative nature of the energy 
and energy fluxes are not vastly different, however, their ampli- 
tudes arc, with the total energy density in the convective region of 
the impermeable case 2-3 orders of magnitude larger than in the 
model bounded by a stable layer. A similar difference in magni- 
tude is seen in the horizontal flux, while the vertical flux is closer 
in the two models. 
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Horizontal Mode 

Figure 5. Kinetic energy spectra for ttie model witli an imper- 
meable lower boundciry (dotted line) and for the model with a 
radiative region at the bottom boundary (solid line). The en- 
ergy is greater at all horizontal wavenumbers in the impermeable 
model. In addition, no dissipation region is seen in the imperme- 
able model, indicating buildup of energy at the smallest scales in 
that case. 
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Figure 6. A snapshot of gravity waves in the radiative interior in 
a completely linear model. Vorticity (/-y^) is shown. This pattern 
is similar to that calculated from ray-tracing. 
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Figure 7. Power spectra (^p-) of gravity waves at the base of the 
convection zone (a) and deeper into the radiative interior (b) / (c) 
for our partially linear model. White represents large energy, blue 
represents low energy. Peak energy in (a) is three orders of mag- 
nitude larger than in (b) / (c) . Energy is broadband in frequency 
and horizontal wavenumber near the base of the convection zone. 
Moving deeper into the radiative interior ridges from standing 
g-modes are seen at high frequencies, while low frequencies arc 
occupied by propagating waves from the overlying convection. 
These spectra were taken shortly after the simulation was begun. 
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Figure 8. Snapshots of our fully nonlinear model. In (a) we show 
the vorticity (/-y^) in the full computational domain. In (b) we 
show the temperature perturbation (white represents positive- 
hotter than reference state temperature, while black represents 
negative-colder than reference state temperature) in the radia- 
tive region only to emphasize the gravity waves. It is seen there 
that the amplitude of the waves falls off quickly and that while 
many scales of g-modes are produced at the interface, only long 
wavelength waves survive in the center. This gravity wave pattern 
is significantly different than that shown in figure 6. 
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Figure 9. Energy spectra, similar to figure 8 but for our fully non- 
linear model and after the convection has reached steady state. 
Within the radiation zone (figures b and c) the peak energy is 
three orders of magnitude larger than in the linear case; the non- 
linear case transfers energy more efficiently to the radiative re- 
gion. However, the high frequency standing g-modes are gone. 
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Figure 10. Temperature perturbation wittiin the radiation zone 
at two different times for liorizontal modes m=20 and 100 arc 
sfiown. The dotted line represents an earlier time, while the solid 
line represents a later time, illustrating the inward movement of 
wave packets and group velocity. The decrease in amplitudes of 
the wave packets shows the radiative damping while moving ra- 
dially inward. 
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Figure 12. Time averaged, real space, kinetic energy density 
/_2m \ function of radius within the radiative interior. A 

^ cms-' ' 

rapid decrease in energy moving radially inward is seen immedi- 
ately, but with some rise in energy closer to the center. 
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Figure 13. Kinetic Energy density as a function of frequency, 
calculated as p{vg solid line and as p{ — )^v^ (estimate from 

the linear dispersion relation) dotted line. Near the convection 
zone the linear dispersion relation significantly underestimates 
the energy density for both modes shown. However, deeper in the 
interior it is a better approximation. 



22 Tamara M. Rogers and Gary A. Glatzmaier 




20 40 60 80 100 120 
Frequency (/-/^Hz) 

Figure 14. Wave energy flux as a function of frequency for m=l 
solid line, ni=10 dotted line and m=20 dashed line, for three 
different radii within the radiation zone. Wave energy flux is dis- 
tributed rather evenly in frequency, particularly for m=l, unlike 
analytic estimates. 



